function [IC_opt,FT_opt] = LyapunovOrbitDetermination_L1_many(u)

% OpenCRTBP_u(u);
hold on
options=odeset('RelTol',1e-14,'AbsTol',1e-22);

[L1,L2,L3,L4,L5] = librationPoints(u);

deltaR = (L2(1)-L1(1))/50;

Ro = L1(1)+deltaR;

InitialCondition_Guess = [Ro;0;0;0;0;0];

FT_Guess = 1.5;

Xo = [-0.0005;FT_Guess] % vy should be negative.

options_fmincon = optimoptions('fmincon', 'MaxFunctionEvaluations', 10000, ...
    'MaxIterations', 50000,'ConstraintTolerance',1e-9,'Algorithm','interior-point');



[Xopt,fval,flag] = fmincon(@CostFunction,...
    Xo,...     % Initial guess
    [], [],... % Linear inequality constraint
    [], [],... % Linear equality constraint
    [], [],... % Boundary constraint
    @(X)DynamicalConst(X,u,InitialCondition_Guess,L1,options),... % Nonlinear Constraint and options.
    options_fmincon);   % Options for fmincon

% Check Initial Condition

IC_opt_small = [Ro;0;0;0;Xopt(1);0];
FT_opt_small = Xopt(2);
% [~,S_opt] = ode113(@(t,S)CR3BP_n(t,S,u),[0,FT_opt_small],IC_opt_small,options);
% S_opt = S_opt';
% plot3(S_opt(1,:),S_opt(2,:),S_opt(3,:),'r','linewidth',3)

% Make Bigger orbit.
IC_temp = IC_opt_small;
FT_temp = FT_opt_small;

IC_opt = [];
FT_opt = [];

while (flag == 1)
    InitialCondition_Guess = [IC_temp(1)+deltaR/5;IC_temp(2:6)];
    Xo(1) = Xopt(1);
    Xo(2) = Xopt(2)
    
    [Xopt,~,flag] = fmincon(@CostFunction,...
    Xo,...     % Initial guess
    [], [],... % Linear inequality constraint
    [], [],... % Linear equality constraint
    [], [],... % Boundary constraint
    @(X)DynamicalConst(X,u,InitialCondition_Guess,L1,options),... % Nonlinear Constraint and options.
    options_fmincon);   % Options for fmincon
    
    if flag ~= 1
        IC_opt_large = [InitialCondition_Guess(1)-deltaR/5,0,0,0,Xo(1),0];
        FT_opt_large = Xo(2);
    end

    

    IC_temp = [InitialCondition_Guess(1);0;0;0;Xopt(1);0];
    FT_temp = Xopt(2);
    
    IC_opt = [IC_opt,IC_temp];
    FT_opt = [FT_opt,FT_temp];
    
%     [~,S_opt] = ode113(@(t,S)CR3BP_n(t,S,u),[0,FT_temp],IC_temp,options);
%     S_opt = S_opt';
%     plot3(S_opt(1,:),S_opt(2,:),S_opt(3,:),'r','linewidth',3)
    
end





%%
    function J = CostFunction(X)
        J = 0;
    end

%%
    function [c,ceq] = DynamicalConst(X,u,InitialCondition_Guess,L1,options)
        vy = X(1);
        FT = X(2);
        IC = [InitialCondition_Guess(1:4);vy;InitialCondition_Guess(6)];
        [~,S] = ode113(@(t,S)CR3BP_n(t,S,u),[0,FT],IC,options);
        S = S';
        plot3(S(1,:),S(2,:),S(3,:))

        c = [S(1,end)-L1(1),1-FT,FT-2];
        ceq = [S(2,end),S(4,end),S(6,end)];
    end

end